Anomalous Scaling and Solitary Waves in Systems with Non-Linear Diffusion 
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We study a non-linear convective-diffusive equation, local in space and time, which has its back- 
ground in the dynamics of the thickness of a wetting film. The presence of a non-linear diffusion 
predicts the existence of fronts as well as shock fronts. Despite the absence of memory effects, 
solutions in the case of pure non-linear diffusion exhibit an anomalous sub-diffusive scaling. Due to 
a balance between non-linear diffusion and convection we, in particular, show that solitary waves 
appear. For large times they merge into a single solitary wave exhibiting a topological stability. 
Even though our results concern a specific equation, numerical simulations supports the view that 
anomalous diffusion and the solitary waves disclosed will be general features in such non-linear 
convective-diffusive dynamics. 



PACS numbers: 47.20.Ft,47.56.+r,47.54.-r,89.75.Fb 

Anomalous diffusion is a topic of great current inter- 
est (see e.g. Refs. [H-Q)- ^ is commonly believed that 
the phenomenon has its roots in correlated structures 
and/or disorder in the medium in which the diffusion oc- 
curs which, in turn, induce memory effects. However, it 
is the aim of the present work to demonstrate that a non- 
linearity in the diffusion process, local in space and time, 
is sufficient to induce anomalous scaling. When, in addi- 
tion, convection enters the non-linear diffusion problem 
we consider, we find, in particular, that the same non- 
linearities that cause the anomalous diffusion, also lead 
to the occurrence of solitary waves exhibiting a topolog- 
ical stability. The topological stability manifests itself 
through the solitary waves having their speed given by 
their average amplitude. 

Even though the non-linear diffusion-convection equa- 
tion that we study has its origins in the dynamics of 
wetting films, related equations appear in other context 
like non-linear heat conduction Wetting films have, 

in their own right, received much attention recently due 
to their relevance for flows in porous media (see e.g. 
Refs.@) and soil physics and hydrology @- Shock- wave 
features are in general expected in systems with non- 
linear convection and diffusion Q . In elastic media such 
features has recently been observed in the laboratory [t| . 
Due to the generic non-linearities in the equations of mo- 
tion an exact mathematical analysis is often not possi- 
ble and numerical simulations are then expected to be 
mandatory. The lack of an exact mathematical analysis 
may, however, prevent one from a deeper insight into the 
physical consequences of the non-linear dynamics. 

In order to be specific, we consider a capillary pore 
flow with a wetting film flowing in a closed wedge with 
a length L as illustrated in Fig[T] The geometry chosen 
appears not be essential but enables us to derive specific 
and analytical predictions. We denote the position along 








Figure 1: Flow of a wetting film in a wedge. The shaded 
area denotes a cross-section area a(x, i) of the film which, in 
general, various along the wedge. 

the wedge with x such that < x < L, and where, e.g., 
the endpoints could be connected to other pores. For 
reasons of simplicity we will consider an infinite half-line 
and hence make L arbitrarily large. The cross sectional 
area of the film at x is a = a(x,t), where t > is the 
time variable. The total and fixed cross sectional area 
of the wedge is At so that the remaining cross section 
is 5A(x, t) = At — a. Apart from the wetting film, the 
closed wedge is filled with a non-wetting fluid. There 
is a pressure drop P = P{x,t) along the wedge in the 
non-wetting fluid and a pressure drop p = p(x, t) in the 
wetting film. The non-wetting fluid flows with a volume 
flow rate Q = Q{x, t) and the wetting fluid with a volume 
flow rate q = q(x, t) so that the total Qt = Q + q is 
constant. Both fluids are assumed to be incompressible. 
The non-wetting fluid has a viscosity M and the wetting 
fluid a viscosity \i. 

Conservation of mass dictates that the area a is related 
to the flow rate through the equation 
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We consider a flow in the films governed by the Darcy 
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equation [lOj], i.e. 



k{a) dp 
[i dx 
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Here k(a) = a 2 /8ir is the effective and a-dependent film 
permeability of the film, assuming a Hagen-Poiseuille 
flow [l(J. The Young relation now gives the connection 
between pressure difference across the interface between 
wetting and non-wetting fluids and the curvature of the 
interface r at capillary equilibrium [To| . i.e. P—p = 7/r, 
where 7 is the surface tension. If the wedge has an open- 
ing angle cj) and the wetting angle is zero, then a geo- 
metrical consideration leads to a = r 2 (cot (cj)/2) — </>/2). 
Combining this relation with the Young relation leads to 
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where we have defined the effective surface tension 7 = 
js/a/r. If we now eliminate p and q between Eqs. (jTJ) 
we find our basic equation of motion 
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We now make the assumption that the pressure in the 
non-wetting fluid, P, is unaffected by the area a of the 
film. This is true if 5 A » a. The Darcy equation 
© then reads MQ/K = -dP/dx, where K = A^/8tt is 
the effective permeability of an assumed Hagen-Poiseuille 
bulk flow. Due to conservation of mass of the non- wetting 
fluid we have dQ/dx — 0, which, after a straightforward 
integration, leads to 



8P(x,t) 
dx 



8MirQ 
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where the pressure drop at x — is neglected, i.e., 
dP(0,t)/dx = 0, and Q ~ Q > is a characteristic 
volume flow rate. All time dependencies enter through 
the non-wetting pressure P at x = 0. If there is no flow 
of non-wetting fluid in the pore, Q = 0, the pressure P 
is constant along it. 

The two fundamental equations in the present paper 
for the evolution of the film in the wedge are Eqs. (0| 
and l[5J). We non-dimensionalize and combine them by 
setting a — aa, x — x£ and t = tt, where a,x,f are 
some typical area, length and time-scales. The combined 
equation is the diffusion- convection equation 



d ( /—da\ da da 

d-A w ^di)- va ^ = d^ 
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The effective diffusivity D — Wy/a is counter intuitive 
when interpreting the diffusion process literary: The dif- 
fusion coefficient increases with increasing concentration 
contrary to conventional diffusion where higher concen- 
tration means a more crowded system and, hence, slower 
diffusion. The non-linear convective term vada/d^ in 
Eq. ([6]) is well known from surface wave propagation in 
shallow water (see e.g. Ref. [8(). It corresponds to a wave 
moving towards increasing £ with a speed proportional to 
a. In the context of shallow-water theory, it is the term 
that produces breaking waves: An initially smooth wave 
will eventually produce a singularity at the leading side. 

When there is no non-wetting fluid flow, only the dif- 
fusive term is present on the left hand side of Eq. (El). 
Eq. ([6]) without convection was considered Rcf.|ll( as 
a model for experimental results of imbibition processes. 
The experimental presence of an imbibition front was the- 
oretically inferred using approximative solutions of Eq. 
(EJ) [5]. Here we point out the existence of such fronts 
as well as shock fronts as exact consequences of Eq. (EJ) . 
It is also a remarkable consequence of our analysis that 
anomalous scaling laws of the form £ 2 ~ V with 7 ^ 1 
can be obtained without invoking fractional derivatives 
[!| in diffusion-like equations (see e.g. Refs. @) or by re- 
ferring to stochastic Mori-Lee equations [l2| and memory 
effects (see e.g. Refs. [HI (HI)- If the effective diffusivity 
D is an increasing (or decreasing) function of a we find 
sub-diffusion with 7 < 1 (or super-diffusion with 7 > 1). 
A related observation of the appearance of anomalous 
diffusion in a class of non-linear Fokker-Planck equations 
can be found in Ref. [l5[ • A renormalization group analy- 
sis of the so called Barenblatt's non-linear diffusion equa- 
tion without convection, which describes the filtration of 
a certain compressible fluid through an elastic porous 
medium, also leads to anomalous diffusion (T^]. The dif- 
fusive non-linearity has then, however, a completely dif- 
ferent origin and form as compared to the non-linearities 
considered in the present paper. 

The non-linear diffusion equation obtained from 
Eq. (EJ) with v = can now be transformed from a par- 
tial differential equation to an ordinary one in terms of a 
scaling function h such that a(£,r) = /i(C//( T )) /f( T )- 
This choice of scaling of the solution is not unique but 
is, by construction, such that a load E defined as E = 

Jq^ x d^a^, r) is conserved in time provided it is finite. 
We now obtain 



as well as 



ho = Vh — — h cyh 
ay 



(7) 
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where w = a 1 / 2 f7/(167r/xx 2 ) and v = 2MafQM/(A? r xfi) 
can be regarded to be the relevant physical parameters in 
the system. The diffusivity in Eq. (E| is D = Wy/a and 
there is also a convective flow velocity va both of which 
therefore are now non-trivial functions of a. 



in terms of two constants, ho and c, of integration. Here 
we have defined the natural Boltzmann scaling variable 
y = £//( T ) We notice that the scaling form of 

the solution and Eq. (© leads to anomalous scaling since 
h ~ constant corresponds to £ 2 ~ r 4 / 5 , i.e. sub-diffusion, 
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Figure 2: Summarizing, with h(0) = 3, the three branches of 
the integral curves Eq. (|10[) with fe = — 1 (curve (i) with a 
shock front), and k = 1 (curve (iii) with a front). In the case 
k = (curve (ii)) we make the replacement \ho\ — > h(0) 5 ^. 



and that solutions of the scaling form a(£, r) = h (£/<7(t)) 
lead to the conventional scaling g(r) ~ r 1 / 2 . The con- 
stant of integration ho in Eq. ||2J) can have either sign. 
The constant c > can be determined from a condition 
on da(^,T)/dr at r = or by the load E. As long as 
ho ^ we can now write the solution a(£, r) in the form 



(\h \c2)V*f(t) 



(9) 



where we have defined y = C c /(I' 1 o|c 2 ) 1 ^ 5 /(t) and where 
h = h(y) now obeys the equation 
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with k = ho/\ho\ = ±1. The scale transformations 
£ -» A£;r -» A 2 t;c -» c/A 2 and |/i | — > \h \/X leaves y 
and h invariant and hence also all solutions a(£, r) in the 
form as given by Eq. (O . Let us now assume that the film 
is concentrated initially around the origin, £ = 0. Hence, 
da/dy < 0, and we must have that k = — 1. Eq. (fTU|) then 
shows that the solution h flows towards ft(i/ max ) = at 
some finite value y = y max with an infinite derivative, 
i.e. we have a shock front. Even though we do not, at 
present, have an analytical result for y max given /i(0), one 
can easily calculate j/ max numerically. We then find that 
?/ max ~ \ogh(0), at least for large values of h(0). Since 
2/m a x is finite the load E will also be finite, which then 
leads to a possible determination of c in terms of E. If, 
on the other hand, k = 1 we find, using Eq. (flU)) . that 
h will have a maximum and then flows to zero asymp- 
totically according to h ~ 1/y. In this case the load 
-E is not well defined and the constant of integration c 
can then be determined using initial data as mentioned 
above. Finally, we have to consider the special case with 
fe = in Eq. (fTUl) which then can be solved explicitly. 



In the scaling variable y one then makes the replacement 
\ho\ — > h(0) 5 / 4 and we can then write 
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where «o = ce(0, 0). In this case the constant of integra- 
tion c can, e.g., be determined by the conserved load E. 
In passing, we notice that the exact solution Eq. (fTTj) cor- 
responds to the known Zeldovich-Raizer solution in the 
context of non-linear heat conduction 0, We show 
the various scaling solutions h(y) of Eq. (fTU|) graphically 
in Fig. [5] for a specific initial value. As we have shown 
above the solutions with fe = — 1 and k — have then 
a finite range. All solutions exhibits anomalous scaling 
due to Eq. ||SJ}. 

We will now turn to the existence of solitary wave prop- 
agation in the pore. In the following we solve the film 
equation Eq. ([5]) under the assumption that the system 
is infinitely long. Hence, we make a shift of the origin 
x = and let L — > oo so that £ can take any value on 
(— oo,+oo). Intuitively, when, in a non-linear system, 
the convective term is non-zero, there may be propagat- 
ing solutions where the convective term and the diffusive 
terms balance each other as is the case for the Burg- 
ers equation [8|. The dynamical equation we consider, 
i.e. Eq.©, actually has the same convective term as the 
Burgers equation but describes the diffusion in a non- 
linear manner. In our case, the primary physical reason 
for the balance alluded to above therefore is the effective 
diffusion constant which is proportional to y/a and hence 
increases the diffusion rate with increasing a. 

By defining the natural variable z — £ — Ut with 
U > , then, after an integration, Eq. ^ is transformed 
into 
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where K is an integration constant. We now rewrite this 
equation in the form 



da 



(a — a + )(a — a_) dz 2w 
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where we notice that = 2K/v and that the velocity 

U of the solitary wave is given by 



a + + a. 



2U 

v 



i.e. 
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We now search for solitary like solutions a s = a to 
Eq. (|13p whose variations are localized in space, i.e., 
dots/dS, — > as |£| becomes large or, equivalently in 
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Figure 3: FEM solution of Eq. © with w = 1.0, v = 
2.0 with an initial two-step configuration a(£, r = 0) = 
[tt/2 - arctan(50£ + 25)] + [tt/2 - arctan(50£ - 25/2)]/2. The 
subsequent evolution of this profile is shown at intervals of 
At = 0.025. The velocity of the profiles at later times, corre- 
sponding to the solitary wave Eq. 1)161) . fulfills Eq. 1)1411 in the 
main text. 



terms of the variable z, as r — > ±oo. In these limits, 
the integration constant K defined in Eq. (|12[) becomes 
K = Ua± — va±/2, where we chose a+ > a_ > 
for the two asymptotic areas of the film and hence 
U > 0. Since a± are real we conclude from Eq. (JTSJ) 
that < K = va + a-/2 < U 2 /2v. Eq. |[T2J) can now be 
integrated in a straightforward manner, i.e., 




— (a+ - a-) (z- 2z ) , 
2w 



(16) 



where zo constitutes the initial condition. We note that 
a s — > a + as t — t oo and a s — > a_ as r — > — oo as long as 
the initial conditions are chosen such that a+ > a s (zo) > 
a_. The left-hand side of Eq. (TTB1 now describes a a s - 
profile that is translated in time. With a given initial 
data a s (zo) the asymptotic values a± are finite. It is a 
straightforward calculation to show that the solution a s 
of the form Eq. (jTfi")) is stable for small deformations at 
large values of |£|. This follows from the observation that 
any solution a = a s + 8 of Eq. (j6|) leads to a diffusion 
equation for 5 keeping terms to linear order in 6 for large 
values of |£|. 

We show in Fig. [3] a finite-element (FEM) solution to 
Eq. (J6j) for an initial arctan-profile with two steps. The 
figure shows the initial profile of a versus £ and subse- 
quent profiles at regular time intervals. We observe that 
the initial profile progressively deforms into a stable pro- 
file that moves at constant speed. This is the solitary 
wave a s as given by Eq. (flo) . Indeed, the initial step to 
the left in Fig. [3] moves with larger velocity than that to 
the right, eventually overtaking and absorbing it. This 
is clear from Eq. (|14[) , which shows that the velocity of a 
solitary wave is proportional to the average a to the left 



Figure 4: FEM solution of Eq. © with w = 1.0 and v = 2.0. 
The initial configuration is a(£, r = 0) = (tt/2) — arctan(50a; — 
25). The subsequent evolution of the profile is shown at inter- 
vals of At = 0.05. The velocity of the profiles, corresponding 
to the solitary wave Eq. ()16p . fulfills Eq. (|14[) in the main text. 



and to the right of it. This is generally true: As long as 
a to the far left is larger than the a to the far right, all 
deformations will be absorbed into a single solitary wave 
moving with the velocity given by Eq. (|14[) . In other 
words, the solitary wave exhibits a topological stability. 
In Fig. 2] we show how an initial single-step configuration 
rapidly evolves into the stable solitary wave a s . 

In conclusion we have shown that the combined ef- 
fects of non-linear convection and diffusion for the dy- 
namics of a wetting film can lead to a rather intrigu- 
ing physics which, in our case, can be analyzed exactly. 
The local dynamical equation used describes a balance 
between viscous and capillary forces at the interface be- 
tween the wetting film and a bulk liquid. When the bulk 
flow rate drops to zero along a pore, perturbations of the 
film spread diffusively and this non-linear diffusivity may, 
as we have seen, lead to the appearance of shock fronts 
and, remarkably, to anomalous diffusion. In the pres- 
ence of non-linear convection in addition to non-linear 
diffusion, we have seen the appearance of solitary waves. 
Numerical analysis has shown that the presence of sev- 
eral wave fronts eventually merges into a single solitary 
wave, the dynamical origin of which we have outlined. 
Concerning potential applications we, e.g., notice that 
drainage processes, where a less wetting fluid displaces 
a more wetting fluid in a porous medium, appears to 
be fairly well understood today. This is, however, not 
the case in the opposite situation which we describe, i.e. 
imbibition, where a more wetting fluid displaces a less 
wetting fluid. We have also observed that the dynami- 
cal equation considered also enters in the context of non- 
linear heat conduction where solitary waves therefore also 
may appear if a suitable form of convection is present. 
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